Transferable relativistic Dirac-Slater pseudopotentials 
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We present a method for constructing a scalar-relativistic pseudopotential which provides exact 
agreement with relativistic Dirac-Slater all-electron eigenvalues at the reference configuration. All- 
electron wave functions are self-consistently computed in the valence region at the exact all-electron 
scalar relativistic eigenvalues. This method improves transferability of the resulting pseudopotential 
and presents a better starting point for the designed non-local pseudopotential approach [Phys. Rev. 
B 59, 12471 (1999)]. We present calculations for the gold atom as an example of the new approach. 
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Ab initio density functional theory |l],|2| (DFT) calcu- 
lations have been used extensively over the past thirty 
years to study metals, insulators and semiconductors. 
The plane- wave pseudopotential method has been widely 
used because it is both accurate and fast, enabling calcu- 
lations on large systems inaccessible with other methods. 
The pseudopotentials in these calculations mimic the ef- 
fect of the nuclei and core electrons on the valence elec- 
trons, dramatically reducing the computational cost of 
the calculations. However, poorly constructed pseudopo- 
tentials will lead to inaccuracies in the solid-state results. 
The most important indicator of the quality of the pseu- 
dopotential is transferability, the ability of the pseudopo- 
tential to reproduce the all-electron valence eigenvalues of 
various atomic configurations. In addition, pseudopoten- 
tials requiring a smaller basis set can significantly speed 
up the solid-state calculation. 

Several methods have been developed to accomplish 
these goals. The ab initio pseudopotentials of Hamann, 
Schliiter and Chiang |}| employ a different spherically- 
symmetric potential for each angular momentum, which 
allows enforcement of the norm-conservation condition 
at one atomic configuration, the reference configura- 
tion. Norm conservation guarantees exact agreement be- 
tween the all-electron and pseudopotential eigenvalues 
and wave functions beyond the core radius r c for the 
reference configuration, and it greatly improves trans- 
ferability in other configurations. The semilocal norm- 
conserving pseudopotentials were further developed by 
Kleinman and Bylander Q] who transformed the pseu- 
dopotential into fully separable nonlocal form, dramat- 
ically reducing the memory cost. The optimized pseu- 
dopotential method developed by Rappe et al. || reduces 
the plane-wave cutoff necessary to use norm-conserving 
pseudopotentials. The designed nonlocal (DNL) ap- 
proach of Ramer and Rappe || significantly improves the 
transferability without affecting the convergence proper- 
ties. Goedecker et al. have developed highly trans- 
ferable dual-space multiple-projector pseudopotentials 
which are expressed as sums of Gaussians. The Van- 
derbilt ultrasoft pseudopotential method || further low- 
ers the plane-wave cutoff by discarding norm conserva- 



tion but restoring the correct charge with augmentation 
charge density functions. In this approach, good transfer- 
ability is achieved by the use of several projectors. The 
combination of ultrasoft pseudopotentials with BlochPs 
projector-augmented wave H method can provide ex- 
tremely good convergence and transferability. A feature 
common to all these methods is the preservation of exact 
agreement of the all-electron and valence eigenvalues at 
the reference configuration. 

Many technologically important materials contain ele- 
ments with Z > 54. For such heavy atoms, relativistic 
effects cannot be ignored without severe consequences 
for the accuracy of the calculations. The most impor- 
tant change to the Kohn-Sham Hamiltonian is in the 
kinetic-energy operator. In DFT all-electron atomic cal- 
culations the Dirac-Slater and the Koelling-Harmon |Tc| ] 
approaches are widely used. The Dirac-Slater equation 
includes spin-orbit coupling, producing (for s = 5) a 
pair of spin-dependent (up and down) wave functions and 
eigenvalues for every orbital with I > 0. The Koelling- 
Harmon approach omits the spin-orbit interaction from 
the Hamiltonian but retains all other relativistic kine- 
matic effects. Therefore a single "spin-averaged" wave 
function and eigenvalue is produced for each atomic or- 
bital. This is less accurate than the Dirac-Slater method, 
but the error thus introduced is often acceptably small. 
In addition to kinetic and spin-orbit effects, the exact 
relativistic correction for the local density approxima- 
tion (LDA) exchange functional is also known but has 
a small effect on the valence states. Recently the rela- 
tivistic exchange functional has been derived within the 
generalized gradient approximation JPH . 

The direct relativistic effects on the valence states of 
a heavy atom are small. Relativity strongly affects the 
core states, changing the self-consistent potential seen 
by the valence orbitals. This implies that a pseudopo- 
tential can effectively incorporate relativistic effects in 
non-relativistic solid-state pseudopotential calculations. 
With a Koelling-Harmon reference atomic all-electron 
calculation, the pseudopotential construction is in no 
way different from the non-relativistic case. However if a 
Dirac-Slater all-electron calculation is used as a basis for 
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the the pseudopotential construction, there exists a prob- 
lem of representing the the up and down relativistic wave 
functions and eigenvalues by a single wave function and 
eigenvalue arising from the non-relativistic Schrodingcr 
equation. Following the suggestion of Klcinman, |12j ] 
Bachelet and Schliiter Jl3| proposed the construction of 
separate V"f( r ) and V^° wn (r) for the up and down states. 
They then create an average pseudopotential weighted by 
the different j degeneracies of the I ± \ states 
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C(0 = ^ TO own W + (I + l)Vj(r)]. (1) 

The Bachelet-Schliiter (BS) averaged pseudopotential 
contains all scalar parts of the relativistic pseudopoten- 
tial. This approach gives rise to an error of order a 2 , 
where a is the fine structure constant (a = 13 y 04 )- How- 
ever the eigenvalues of the BS averaged pseudopotential 
will not be equal to the weighted average of the up and 
down all-electron eigenvalues. In addition, the valence 
charge density due to the solutions of the Schrodingcr 
equation for the BS averaged pseudopotential will not 
be the same as the valence charge density of the all- 
electron atom. Therefore, exact agreement between the 
pseudopotential and the all-electron orbital eigenvalues 
and charge density past r c is not preserved even for the 
reference configuration. Typically the error in the eigen- 
values and total energies introduced by the potential av- 
eraging is small, about 1-4 mRy. None of the modern 
pseudopotential methods can remove this error. Thus 
one currently has a choice of either accepting the error 
due to ignoring the spin-orbit coupling in the all-electron 
calculations (Koelling-Harmon) , or using the more cor- 
rect Dirac-Slater equations and accepting the error due 
to pseudopotential averaging. 

We present below a method for constructing a pseu- 
dopotential whose eigenvalues agree exactly with the av- 
eraged relativistic all-electron eigenvalues 
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It is impossible to enforce charge density agreement with 
the relativistic all-electron calculation for all r > r c , since 
each wave function asymptotes to an exponential, and 
the sum of two exponential functions cannot be equal to 
a single exponential function at more than two points. 
Instead, we impose an aggregate norm conservation out- 
side the core radius r c . We require the weighted average 
of the integrals of the squares of the spin-up and spin- 
down wave functions from r c to infinity to be preserved 
by the new potentials and pseudo-wave functions for each 
valence orbital 4> n i 
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Unlike the standard averaging procedure where the 
spin-up and spin-down pseudopotential s are averaged 
and the averaged pseudopotential is then used to obtain 
the non-relativistic wave function and eigenvalues, in our 
procedure the relativistic all-electron eigenvalues and the 
aggregate norm conservation criterion are used to obtain 
the non-relativistic wave functions and potentials at the 
all-electron level. 

The algorithm proceeds as follows. First we obtain the 
relativistic spin-up and spin-down all-electron wave func- 
tions and eigenvalues using the Dirac-Slater formalism. 
The averaged all-electron eigenvalue and norm, ef^ G and 
Q n i, are then computed for each valence orbital. The ini- 
tial guess for the new non-relativistic 4> n i is set to be the 
weighted average of the spin up and spin down Dirac- 
Slater wave functions ip^f and ^f wn , scaled to preserve 
Qni. The valence charge density p V aiencc(f) correspond- 
ing to ip^f and ipn° wn is computed and then added to the 
charge density due to the core orbitals PcorcM to give 
total charge density ptot&\(r). 

The exchange-correlation and Hartree potentials and 
energy densities are then computed and the full potential 
V(r) is obtained. For each orbital the potential V(r) is 
used in an inward solve ]l4[ ] for the new wave function 
<f) n i with the criterion e„; = e^ G - Each orbital is scaled 
so that the norm outside r c is equal to Q n i- For r < r c , 
4> n i is given a smooth nodeless form 



AAVG 



(r c ) 



Cnl 1 - 



(4) 



with the parameter c n i chosen such that 4> n i is normalized 
to unity. This is done for each valence orbital. The new 
valence charge density p va .ience(r) is then computed from 
the 4>ni , completing the cycle. This cycle is repeated until 
the input and output <f> n i are equal. 

Pseudopotentials arc then generated from the con- 
verged 4> n i. The optimized and designed nonlocal pseu- 
dopotential methods are used in this paper to improve 
convergence and transferability. However we empha- 
size that the present method of insuring agreement be- 
tween relativistic all-electron and pseudopotential re- 
sults can be used with any pseudopotential generating 
scheme, such as the ultrasoft construction, since the self- 
consistent solve is done at the all-electron level. 

We have applied our new averaging method to the gold 
atom. Gold has atomic number 79 and therefore must 
be treated relativistically. We have used a neutral 6s 10 
6p 5c? 10 reference configuration. We use the Perdew- 
Zunger [jl6| parametrization of the Ceperley-Alder (ll| 
LDA for the exchange-correlation functional in all calcu- 
lations. 

As can be seen from the results in Table I, the new 
method improves transferability in configurations close 
to reference, with the sum of absolute values of er- 
rors in orbital eigenvalues of neutral configurations re- 
duced by 60%. Even more importantly, rigorously cor- 
rect eigenvalues for the reference configuration enable the 
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designed nonlocal method to improve transferability fur- 
ther. Since the DNL augmentation operator does not 
affect the pseudo-wave functions and eigenvalues of the 
reference configuration, the error given by the pseudopo- 
tcntial averaging method can never be corrected with the 
designed nonlocal formalism. Furthermore, the augmen- 
tation operator has less impact on configurations close 
to reference than on those further away; thus eigenvalue 
errors in the neutral configurations are much harder to 
eliminate for the pseudopotential averaging method than 
for our new all-electron averaging method. While the to- 
tal transferability error for the new all-electron averaging 
based designed nonlocal potential is 40% less than that 
of the BS averaging based designed nonlocal potential, 
the error in neutral configurations is 78% less. 

We have presented in this paper a more accurate 
method for construction of relativistic pscudopotentials 
based on Dirac-Slater all-electron atomic calculations. 
Our procedure enforces the agreement of the pseudopo- 
tential and the all-electron spin averaged eigenvalues 
at the reference configuration, leading to a significant 
improvement in pseudopotential transferability. The 
method has negligible computational cost, is easy to im- 
plement, and applies to both norm-conserving and ultra- 
soft pseudopotential formalisms. 
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tory for Research on the Structure of Matter and the 
Research Foundation at the University of Pennsylvania. 
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Supercomputer Center and the National Center for Su- 
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TABLE I. Configuration testing for the Bachelet-Schliiter 
(BS) Au pseudopotential averaging and our all-electron aver- 
aging methods described in text. Eigenvalues and eigenvalue 
errors are given for all-electron (AE) weighted average eigen- 
values, BS averaged and our all-electron averaged pseudopo- 
tential eigenvalues with and without the designed nonlocal 
(DNL) construction. All energies are in Ry. 





AE 
average 


BS 
averaging 
error 


Present 
averaging 
error 


BS 
with DNL 
error 


Present 
with DNL 
error 


6s 1 
6p° 
5d 1() 


-0.4457 
-0.0696 
-0.5281 


0.0011 
0.0035 
0.0035 


-0.0000 
0.0000 
0.0000 


0.0011 
0.0035 
0.0035 


0.0000 
0.0000 
0.0000 


6s 2 
6p° 
5d 9 


-0.5096 
-0.0968 
-0.6783 


0.0013 
0.0050 
-0.0027 


0.0001 
0.0006 
-0.0064 


0.0018 
0.0050 
0.0000 


0.0010 
0.0007 
-0.0019 


6s° 
6p° 
5d w 


-0.9958 
-0.5118 
-1.1509 


-0.0008 
0.0038 
0.0037 


-0.0013 
-0.0037 
-0.0009 


0.0009 
0.0036 
0.0041 


-0.0009 
-0.0036 
0.0006 


6s 1 
6p° 
5d 9 


-1.0786 
-0.5653 
-1.3211 


0.0010 
0.0059 
-0.0031 


-0.0011 
-0.0024 
-0.0076 


0.0016 
0.0052 
0.0005 


0.0007 
-0.0018 
-0.0009 


6s 2 
6p° 
5d 8 


-1.1611 
-0.6152 
-1.5009 


0.0009 
0.0081 
-0.0099 


-0.0012 
-0.0010 
-0.0143 


0.0025 
0.0064 
-0.00183 


0.0023 
-0.0003 
-0.0022 


6s° 
6p° 
5d 9 


-1.6971 
-1.0960 
-2.0317 


0.0022 
0.0052 
-0.0041 


-0.0026 
-0.0065 
-0.0093 


0.0007 
0.0015 
0.0001 


0.0001 
-0.0055 
-0.0001 


6s 1 
6p° 
5d 8 


-1.7974 
-1.1694 
-2.2275 


-0.0018 
0.0074 
-0.0111 


-0.0028 
-0.0048 
-0.0160 


0.0012 
0.0010 
-0.0021 


0.0015 
-0.0030 
0.0013 


6s 2 
6p° 
5d 7 


-1.8947 
-1.2362 
-2.4296 


-0.0015 
0.0100 
-0.0154 


-0.0026 
-0.0028 
-0.0201 


0.0026 
-0.0010 
0.0002 


0.0033 
-0.0002 
-0.0001 
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